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1.0  INTRODUCTION 


This  report  summarizes  the  research  accomplishments  achieved  under  Grant  AFOSR 
F49620-00- 1-0229,  “Application  of  Large  Eddy  Simulation  to  Cooling  and  Flow 
Problems  in  Aeropropulsion  Systems.”  The  objective  of  the  research  was  to  expand  the 
capabilities  of  LES  technology  to  contribute  to  the  solution  of  urgent  problems  in 
propulsion  systems  and  to  contribute  to  the  physical  understanding  of  such  flows.  The 
long-term  goal  of  the  thrust  is  to  advance  large  eddy  simulation  (LES)  technology  to  the 
point  where  it  is  effective  for  computing  realistic  flows  in  a  rotating,  multistage 
turbomachine.  First,  however,  some  building  block  simulations  need  to  be  carried  out  for 
a  single  blade  element.  The  first  task  in  the  long-term  research  was  to  address  cooling 
flows  in  turbines,  both  the  film  cooling  of  external  blade  surfaces  and  the  complex  flows 
in  internal  passages. 

The  grant  duration  was  a  short  eight  months  so  the  scope  of  the  work  that  could  be 
completed  in  that  period  was  very  limited.  Work  was  initiated  on  the  external  film 
cooling  problem  and  some  interesting  results  were  obtained  from  a  study  on  the  effect  of 
rotation  on  heat  transfer  in  channel  flow.  The  rotating  channel  flow  heat  transfer  results 
represent  one  step  toward  the  realistic  simulation  of  internal  cooling  flows  for  turbine 
components.  Rotation  exerts  an  important  influence  on  flows  in  turbomachines  and 
needs  to  be  taken  into  account  in  the  design  process.  The  remainder  of  this  report  will 
describe  results  obtained  from  the  rotating  channel  study. 

System  rotation  influences  turbulence  in  several  ways:  for  example,  it  may  decrease 
energy  transfer  from  large  to  small  scales  or  reduce  turbulence  dissipation  and  decay  rate 
of  turbulence  energy.  In  rotating  channel  flow,  system  rotation  can  both  stabilize  and 
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destabilize  the  flow.  On  the  unstable  side,  Coriolis  forces  enhance  turbulence  production 
and  increase  the  intensity  of  turbulence,  while  on  the  stable  side,  Coriolis  forces  reduce 
turbulence  production  and  decrease  the  intensity  of  turbulence.  Most  large  eddy 
simulations  (LES)  of  flows  with  heat  transfer  reported  to  date  employed  an 
incompressible  formulation  and  treated  temperature  as  a  passive  scalar.  However, 
compressible  formulations  have  been  employed  in  a  few  recent  works  (Moin  et  al. 
(1991),  Nicoud  (1998),  Wang  and  Fletcher  (1995;  1996),  Dailey  and  Fletcher  (1996; 
1997;  1998),  and  Meng  et  al.  (1999a;  1999b))  where  the  coupling  between  velocity  and 
temperature  fields  were  considered. 

Many  rotational -induced  flow  phenomena  have  been  reported  on  by  Tritton  (1978; 
1985)  and  Hopfinger,  (1989).  Rotating  channel  flows  have  been  investigated- 
experimentally  (Johnston  et  al.,  1972;  Han  et  al.,  1994)  and  numerically  (Kristofferson 
and  Andersson,  1993;  Tafti,  1991)  (  by  direct  and  large-eddy  simulations).  However, 
very  little  information  is  available  for  rotating  channel  flow  with  heat  transfer,  especially 
from  the  LES  or  DNS  community.  The  goal  of  this  facet  of  the  research  was  to 
demonstrate  and  analyze  the  effects  of  system  rotation  on  heat  transfer  in  turbulent 
channel  flow. 

2.0  GOVERNING  EQUATIONS 

The  governing  equations  for  large  eddy  simulation  were  obtained  by  filtering  the 
unsteady  three-dimensional  compressible  Navier-Stokes  equations.  The  filter,  which 
separates  the  effects  of  the  large-scale  and  small-scale  motions,  can  be  written  in  terms  of 
a  convolution  integral  as 


/(x,0  =  J^G(x-|)/(x,04 


(1) 
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where  G  is  some  spatial  filter  with  a  width  on  the  order  of  the  grid  spacing,  and  D  is  the 
flow  domain.  For  compressible  flows,  the  Favre-averaged  filter  is  introduced  which  can 
be  defined  by 

jM  (2) 

P 


The  resulting  equations  were  nondimensionalized  with  respect  to  a  reference  length 
6,  velocity  V,,  density  Pr,  and  viscosity  jXr-  For  the  simulations  in  this  report,  5  is  the  half 
height  of  the  channel,  and  Vr  is  the  bulk  velocity  at  the  inlet.  The  resulting  Favre-filtered 
nondimensional  compressible  forms  of  equations  for  conservation  of  mass,  momentum 
are  given  below  : 


ap  ^  a(pu,)_Q 

0/  dx , 


(3) 


a(pH,)  ^(p«,«,)_  dp  dOy  dX^j 


^(P£)  +  -^|p£  +  p)mJ=- 

dt  dx. 


dxj  dxj  dxj 


-a-n-e 


(4) 


(5) 


In  the  equations  above,  2pey^Q,Uj  represents  the  Coriolis  force,  8^3  is  Levi-Civita's 

alternating  tensor  and  Q  is  the  angular  velocity  of  the  system.  A  traditional  right-handed 
coordinate  system  is  employed  and  the  axis  of  rotation  is  in  the  positive  z  or  X3  direction. 
The  resolved  total  energy  is  £  =  c^f  +\u.u. ,  the  viscous  stress  tensor  is 


ft  =2^ 

'  Re, 


(6) 


where  5ij  is  the  Kronecker  delta  and  the  strain  rate  tensor  and  heat  flux  vectors  are 
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c 

2]^  dxj  ax,.  ^ 

Cf,]i  af 

Re^  Pr  ax^. 


(7) 

(8) 


The  Favre  filtered  equation  ui  state  is  p  =  pRT  .  The  viscosity  and  thermal  conductivity 
were  evaluated  by  using  the  power  law  variation  with  temperature. 


iL 

1^0 


0.7 

k  _ 

^  rj,  ^ 

i^oj 

’  K~ 

i^oj 

0.7 


(9) 


The  Prandtl  number  was  assumed  to  be  constant  at  a  value  of  0.71. 

The  effects  of  the  small  scales  are  present  in  the  above  equations  through  the  SGS 
stress  tensor  in  the  filtered  momentum  equation 


7/,  =  P 


(10) 


The  SGS  heat  flux  and  other  quantities  in  the  filtered  energy  equation  are  given  by 


q,j  =  pc.{Tuj-fu^) 

(11) 

9(pt,^) 

a  =  M,  ' 

dx^. 

(12) 

dw,  _  dw, 

(13) 

d«,  _  du, 

c  =  rr  - i - rr  - L 

dXj  dXj 

(14) 
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For  this  research,  the  first  three  terms  on  the  right-hand  side  of  the  energy  equations  were 
computed,  but  the  remaining  terms  (a,  n,  e)  were  neglected  due  to  their  relative  smaller 
magnitudes. 

3.0  NUMERICAL  METHOD 

In  the  present  research,  a  compressible  finite  volume  formulation  CDailey,  1997; 
Dailey  and  Fletcher,  1998;  Meng  et  al.,  1999b)  was  used  to  solve  the  governing 
equations,  which  also  permits  the  subgrid-scale  turbulent  Prandtl  number  to  be  computed 
dynamically.  The  code  used  Cartesian  hexahedral  control  volumes,  and  solved  for  the 
primitive  variables  (p,  u,  v,  w,  T)  that  were  stored  at  the  cell  centers.  Besides  the  physical 
time  integration,  pseudo  time  iterations  were  performed  to  resolve  nonlinearities  in  the 
algebraic  formulation  and  to  implement  the  low  Mach  number  preconditioning  using  an 
implicit  LU-SGS  scheme.  The  low  Mach  number  preconditioning  was  used  to  enable  the 
compressible  code  to  work  efficiently  at  nearly  incompressible  speeds.  The  solver  is 
second-order  accurate  in  both  space  and  time  (Dailey  and  Fletcher,  1996;  Dailey  1997). 
The  multiblock  code  was  parallelized  using  the  message  passing  interface  (MFI).  The 
computations  were  carried  out  on  an  IBM  SF2  (Minnesota  Supercomputing  Institute) 
using  17  processors. 

The  simulations  were  run  with  a  domain  size  of  2nx2x7iwith  a  grid  that  had  48  x 
64  X  48  control  volumes  in  the  x,  y,  and  z  directions,  respectively.  A  grid  study  (Dailey, 
1997)  has  shown  that  for  low  heating  simulations  without  rotation,  this  grid  size  provided 
accurate  results  compared  to  DNS  and  experimental  data.  The  near  wall  region  was  well- 
resolved.  For  the  high  heating  case,  the  y"^  (based  on  averaged  friction  velocity)  at  the 
first  grid  point  near  the  wall  was  less  than  0.5,  so  no  wall  function  was  used. 
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When  property  variations  are  taken  into  account,  flows  with  heating  or  cooling  do  not 
attain  a  fully  developed  state.  However,  experiments  show  that  far  downstream  of  the 
entry  region,  a  slowly  evolving  quasi-developed  state  exists.  Dailey  (1997)  and  Dailey 
and  Fletcher  (1998)  showed  that  for  laminar  flow  with  uniform  heat  flux  conditions,  use 
of  step  periodic  boundary  conditions  can  match  the  flow  conditions  thai  6Volve  from 
variable  property  developing  flow.  Inspection  of  a  developing  flow  solution  revealed  that 
the  following  conditions  held  approximately  in  the  downstream  region: 


^(p“)_o 

dx 

dx 


=  0 


dw 
dx 
dp 

— -constant 
dx 

dT 

=  constant 


dx 

Thus,  Dailey  (1997)  concluded  that  a  short  section  of  the  downstream  region  could 
be  computed  in  a  "stepwise  periodic"  manner  with  the  following  stream'  /ise  boundary 
conditions. 


pu{0,y)  =  pu{L^,y) 
v(0,y)  =  v(4.y) 
w{0,y)  =  w{L^,y) 

Pp  {0,y)  =  Pp{L,,y) 

r(0.y)  =  7(4.y)-Ar, 


where  Lx  is  the  length  of  the  computation  domain  in  the  streamwise  direction  and  is 


the  periodic  component  of  the  pressure  p{^x,y,z,t)  =  f>x  +  Pp{^x,y,z,t)  where  P  is  the 


average  streamwise  pressure  gradient.  The  temperature  difference  AT  is  computed  from 
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an  energy  balance  utilizing  the  uniform  heat  flux  imposed  and  the  mass  flux.  All 
primitive  variables  were  assumed  to  be  periodic  in  the  spanwise  (z)  direction. 

No-slip  velocity  and  zero  normal  pressure  gradient  boundary  conditions  were 

enforced  at  the  upper  and  lower  walls.  The  isoflux  thermal  wall  boundary  conditions 

♦ 

q 

were  used  for  both  walls.  The  dimensionless  wall  heat  flux,  q^.  = -  .  was  kept 

ref^  ref^  ref 

at  2  X  10'^  for  the  low  heating  case,  and  kept  at  2  x  10’^  and  -2  x  10'^  for  the  high  heating 
and  high  cooling  cases,  respectively. 

3.1  Dynamic  Subgrid-Scale  Stress  Model 
The  compressible  formulation  for  the  dynamic  model  first  proposed  by  Moin  et  al. 
(1991)  was  utilized  in  this  research.  The  dynamic  subgrid  stress  model  uses  information, 
from  two  different  grid  levels  to  determine  the  coefficient  in  the  Smagorinsky  stress 
model  dynamically.  The  dynamic  model  does  not  require  ad  hoc  damping  near  solid 
walls  and  enables  the  turbulent  Prandtl  number  to  also  be  computed  dynamically. 

The  dynamic  model  coefficient  and  the  turbulent  Prandtl  number  were  averaged 
along  the  only  homogeneous  (spanwise)  direction.  Further  details  can  be  found  in  Meng 
and  Pletcher  (2000). 

4.0  RESULTS 

Caculations  were  performed  for  four  different  cases:  no  heat  transfer,  low  heating, 
high  heating,  and  high  cooling.  The  Reynold  number  Re  =  V^25  /v  (based  on  the  channel 

width,  26,  and  reference  velocity,  Vr)  was  5600,  and  the  rotation  number  =|^2|26/V, 
was  0.144.  It  should  be  noted  that  the  results  to  be  presented  here  are  for  a  channel  with 
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negative  (Q  <  0)  angular  velocity  about  the  positive  z-axis.  The  configuration  can  be 
seen  in  Fig.  1. 

It  should  be  noted  that  the  results  to  follow  are  for  channel  flow  rather  than  duct  flow 
so  no  mean  secondary  flow  is  present.  Such  conditions  represent  the  flow  in  the  central 
part  of  a  duct  of  large  (in  the  limit,  infinite)  aspect  ratio.  The  results  from  first  case  were 
compared  with  the  DNS  results  (grid  96  x  97  x  128)  of  Piomelli  and  Liu  (1995)  for  the 
isothermal  rotating  channel  flow  with  Re^  =  i7j,2<5  /v  =  5700  and  =  0.144 .  Since  it 
was  found  from  previous  LES  simulations  that  the  contribution  of  the  SGS  model  to  the 
turbulent  flux  is  small,  it  was  of  interest  to  do  a  simulation  (coarse  grid  DNS)  without  any 
model  using  the  same  grid  resolution  as  that  of  LES  (48  x  64  x  48). 

Figure  2  shows  the  profile  of  streamwise  velocity  versus  y/8.  The  DNS  data 
(Piomelli  and  Liu,  1995)  (replotted  for  negative  angular  velocity)  are  represented  by 
circle  symbols,  the  LES  data  are  represented  by  a  solid  line,  and  the  coarse  grid  DNS  data 
are  represented  by  a  dashed  line.  Good  agreement  has  been  found  between  the  LES  and 
DNS  data.  It  is  interesting  that  the  U  velocity  distribution  is  no  longer  symmetric  about 
the  center  line.  The  maximum  U  velocity  is  shifted  from  the  centerline  toward  the 
stabilized  wall  due  to  the  Coriolis  force.  The  coarse  and  fine  grid  DNS  data  are  almost 
indentical. 

Figures  3, 4,  and  5  show  the  rms  distributions  (normalized  by  the  average  value  of  Uj 
on  two  walls)  for  velocity  fluctuations  in  x,  y,  and  z  directions,  respectively.  The  rms 
fluctuations  are  enhanced  near  the  unstable  side  (y  =  1.0)  but  reduced  near  the  stable  side 
of  the  channel  (y  =  -1.0).  It  is  noticeable  that  the  profile  of  v^s  becomes  a  one-peak 
distribution  instead  of  the  two-peak  distribution  obtained  without  rotation.  This  is  also 
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Fig.  1  Coordinate  system 
for  rotating  channel  flow 


Fig.  3  Streamwise  velocity  rms 
distribution,  isothermal  flow 


y/S 

Fig.  5  Spanwise  velocity  rms 
distribution,  isothermal  flow 


1.5 


Fig.  2  Streamwise  velocity 
distribution,  isothermal  channel  flow 


Fig.  4  Normal  velocity  rms 
distribution,  isothermal  flow 


due  to  the  fact  that  the  Coriolis  force  acts  in  the  y  direction.  These  results  agree 
reasonably  well  with  those  of  Piomelli's,  validating  the  current  LES  formulation.  The 
coarse  grid  DNS  data,  in  general,  slightly  overestimate  velocity  fluctuations  compared 
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with  the  LES  data.  There  is  an  overshoot  for  the  coarse  grid  DNS  in  the  prediction  of 
velocity  fluctuations  near  the  upper  (top)  wall.  The  reason  might  be  that  the  current  grid 
resolution  is  not  fine  enough  for  DNS. 

Figure  6  shows  the  streamwise  velocity  profiles  for  the  high  heating  (HH),  high 
cooling  (HC)  and  low  heating  (LH)  cases.  The  influence  of  rotation  is  reflected  in  the 
asymmetric  distribution  of  velocities  for  all  three  cases.  The  U  velocity  profile  for  HC  is 
above  that  of  LH,  but  the  profile  for  HH  is  below  that  of  LH.  This  is,  in  part,  because  in 
the  figures  the  streamwise  velocity  has  been  normalized  by  the  average  of  Ut  at  the  two 
walls  and  is  significantly  influenced  by  property  variations  in  the  flow.  Previous 
researchers  (Huang  et  al.,  1995;  Dailey,  1997;  Dailey  and  Pletcher,  1998)  have  proposed 
a  semi-local  definition  of  friction  velocity,  in  which  the  local  density  is  used  instead  of 
the  wall  density. 


When  m’  is  used  to  replace  Mx,  the  resulting  semi-lncal  distribution,  seen  in  Fig.  7  shows 

that  the  three  profiles  collapsed  toward  the  low  heating  curve. 

The  mean  temperature  profiles  are  plotted  in  global  coordinates  in  Fig.  8.  The 
friction  temperature  parameters  is  defined  below, 

T,=-^  (15) 

As  was  observed  for  the  velocity  profiles,  the  high  heating  and  high  cooling  temperature 
profiles  nondimensionalized  by  wall  values  depart  from  the  low  heating  results. 


-TVT 
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Fig.  7  Profile  of  semilocal  u 
velocity. 


The  root-mean-square  of  the  velocity  fluctuations  with  respect  to  Rf'y.iolds  ensemble 
averages  are  plotted  in  global  coordinates  in  Figs.  9  and  10.  Again,  larger  velocity 
fluctuations  were  found  near  the  unstable  side  as  a  consequence  of  rotation,  and  the  high 
heating  and  cooling  values  varied  significantly  from  the  low  heating  result.  As  before, 
the  values  are  normalized  by  the  average  of  for  the  two  walls. 

Figures  11,  12,  and  13  show  the  viscous,  resolved,  and  modeled  SGS  shear  stress 
distributions  normalized  by  the  average  of  the  two  wall  shear  stresses.  The  results 
indicate  that  rotation  increased  the  shear  stress  near  the  unstable  side  but  suppressed  the 
shear  stress  near  the  stable  side. 


12 


Fig.  10  Distribution  of  Vnns 


Fig.  12  Distribution  of  resolved 
stress 


Fig.  1 1  Distribution  of  viscous 
stress 


Fig.  13  Distribution  of  SGS 
stress 


It  can  be  observed  that  the  shear  stress  on  the  unstabilized  wall  is  approximately  twice  as 
large  as  on  the  stabilized  wall  and  that  the  resolved  turbulent  stress  follows  a  similar 
trend. 

Figures  14,  15,  and  16  show  the  heat  conduction,  resolved  turbulent  heat  flux  and 
modeled  SGS  heat  flux  distributions  normalized  by  the  wall  heat  flux.  The  same  trends 
can  be  observed  as  for  the  shear  stress  distributions,  suggesting  that  heat  transfer  was 
increased  near  the  unstable  side  but  decreased  near  the  stable  side  because  of  rotation. 
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Fig.  14  Distribution  of 
molecular  heat  conduction 


y/5 


Fig.  16  Distribution  of  SGS 
heat  flux 


Fig.  15  Distribution  of  resolved 
heat  flux 


Fig.  17  Contours  of  wall  normal 
vorticity,  high  heating 


Figure  17  shows  the  wall  normal  instantaneous  y-component  of  vorticity  at  planes 
located  at  y/5  =  ±0.99  for  the  HH  case.  Near  the  trailing  (top)  wall  small  scale  turbulent 
structures  are  evident  but  few  can  be  seen  near  the  leading  wall.  Similar  results  were 
observed  for  all  the  rotation  cases  considered  although  the  magnitudes  of  the  vorticity 


was  generally  smaller  for  the  low  heating  case. 

Friction  and  heat  transfer  parameters  are  given  in  Table  1.  The  tabulated 
AT  =  is  the  temperature  difference  across  the  channel  and  Tb  is  the  bulk 

temperature.  The  Nusselt  numbers  are  based  on  the  wall  heat  flux,  channel  hydraulic 
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Table  1.  Summary  of  Wall  Parameters,  Rotating  Flow 


Case 

AT 

T, 

^  fjop 

Low 

Heating 

0.0696 

0.0468 

-0.0321 

64.27 

19.36 

0.00557 

0.00238 

High 

Heating 

0.0889 

O.GoOl 

-0.8038 

83.04 

15.38 

0.00617 

0.00288 

High 

Cooling 

0.0760 

0.0384 

0.3645 

51.81 

19.73 

0.00519 

0.00248 

diameter,  the  bulk  thermal  conductivity,  and  the  difference  between  the  wall  and  bulk 
temperatures.  The  friction  coefficient  is  defined  as 


The  ratio  of  friction  velocity  at  the  top  wall  to  that  at  the  bottom  wall  is  approximately' 
1.8,  1.5,  and  1.1  for  high  cooling,  low  heating,  and  high  heating,  respectively.  The  low 
heating  case  had  a  relative  temperature  difference  of  3.2%,  but  high  heating  and  high 
cooling  cases  exhibited  much  larger  temperature  differences,  namely,  80%  and  34%, 
respectively.  The  Nusseh  number  at  the  top  wall  was  found  to  be  2-5  times  larger  than 
on  the  lower  wall,  whereas  the  friction  coefficients  at  the  top  wall  were  about  twice  as 
large  as  those  at  the  bottom  wall.  Such  a  large  increase  in  Nusselt  number  due  to  rotation 
for  the  heating  case  was  not  expected  and  the  cause  is  currently  under  investigation. 
Density  is  a  factor  in  the  Coriolis  force  and  significant  density  variations  across  the  flow 
may  well  be  responsible  for  the  sensitivity  of  the  Nusselt  number  to  heating  level. 
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5.0  CONCLUSIONS 

The  LES  results  (mean  velocity  and  velocity  fluctuations)  agreed  well  with 
Piomelli's  DNS  data  for  the  incompressible  isothermal  rotating  channel  flow  at  a  rotation 
number  of  0.144.  Turbulent  channel  flow  with  low  heating  and  with  heating  and  cooling 
rates  of  magnitudes  large  enough  to  cause  significant  variation  in  the  temperature- 
dependent  fluid  properties  were  simulated.  All  simulations  were  performed  under  the 
influence  of  spanwise  system  rotation.  The  distribution  of  mean  velocity  and  mean 
temperature  were  influenced  by  system  rotation,  creating  asymmetric  profiles.  The  shift 
of  peak  values  was  toward  the  stable  side  of  the  channel.  This  shift  would  likely  be 
altered  in  ducts  of  finite  aspect  ratios  where  mean  secondary  flows  would  influence  the 
distributions.  In  general,  the  system  rotation  was  found  to  suppress  turbulent  velocity- 
fluctuations  and  shear  stresses  near  the  stable  side  of  the  channel,  but  enhance  them  near 
the  unstable  side.  Accordingly,  turbulent  temperature  fluctuations  and  turbulent  heat  flux 
are  decreased  near  the  stable  side  of  the  channel,  but  increased  near  the  unstable  side  of 
the  channel.  The  streamwise  and  spanwise  mis  velocity  fluctuations  were  seen  to  differ 
by  about  a  factor  of  two  near  the  two  walls.  Ihe  ratio  of  wall  shear  stress  and  heat  flux 
on  the  two  sides  ranged  between  a  factor  of  2  and  3.  The  ratio  of  Nusselt  numbers  at  the 
two  walls  was  as  large  as  a  factor  of  5  as  can  be  seen  from  Table  1.  Using  semi-local 
wall  coordinates  for  velocity  and  temperature  distributions  tended  to  collapse  curves  of 
high  heating  and  high  cooling  into  that  of  low  heating.  It  is  interesting,  too,  that  the 
distribution  of  Vms  is  changed  significantly  (from  a  two-peak  profile  to  a  one-peak 
profile)  due  to  the  Coriolis  force  acting  in  the  y  direction. 
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